###########################################################
#     Prj: WGCNA all in one
#     Assignment: datExpr genarator and power selection
#     Author: Shawn Wang
#     Date: Mar 24,2021
###########################################################

# options and package -----------------------------------------------------
load("~/02.MyScript/OneStepWGCNA/04.OneInAll/functions.Rdata")
suppressMessages(library(DESeq2))
suppressMessages(library(ggplot2))
suppressMessages(library(stringr))
suppressMessages(library(ape))
suppressMessages(library(reshape2))
suppressMessages(library(edgeR))
suppressMessages(library(ggprism))
suppressMessages(library(patchwork))
suppressMessages(library(tidyverse))
suppressMessages(library(WGCNA))
suppressMessages(library(getopt))
suppressMessages(library(ggprism))
suppressMessages(library(ggtree))
allowWGCNAThreads()
options(stringsAsFactors = F)
# args --------------------------------------------------------------------
command=matrix(c(
  'help', 'h', 0, 'logic', 'help information',
  'rawdata', 'r', 1, 'character', 'inputfile: readcount matrix, geneID in row, sample names in column',
  'RcCutoff', 'c', '2', 'integer', 'Noise remove: based on WGCNA FAQ, background noise should be removed, set the readcount value cutoff.(default = 6)',
  'samplePerc', 'p', '2', 'double', 'Noise remove: At least how many samples have readcount value greater than cutoff,(Range 0-1, default = 0.8)',
  'datatype', 'd', 2, 'character', 'datatype: project name for output file names.(default = system date)',
  'method', 'm', 1, 'character', 'working directory, you have type -w $PWD or -w `pwd`',
  'RemainGeneNum', 'g', 1, 'integer', 'gene num cutoff',
  'cutmethod', 't', 1, 'character', 'MAD or Var',
  'rscut', 's', 1, 'double' , "Rsquare cutoff"
),byrow = T, ncol = 5)
args = getopt(command)
## default valuels==========
if (!is.null(args$help)) {
  cat(paste(getopt(command, usage = T), "\n"))
  #  q(status=1)
}

rawdata = args$rawdata

## default value
if (is.null(args$RcCutoff)){
  args$RcCutoff = 6
}
RcCutoff = args$RcCutoff

if (is.null(args$samplePerc)){
  args$samplePerc = 0.8
}
samplePerc  =args$samplePerc 

if (is.null(args$datatype)){
  args$datatype = "count"
}
datatype =  args$datatype

if (is.null(args$RemainGeneNum)){
  args$RemainGeneNum <- 12000
}
RemainGeneNum = args$RemainGeneNum
if (is.null(args$method)){
  args$method <- 12000
}
method = args$method

if (is.null(args$cutmethod)){
  args$cutmethod <- "MAD"
}
cutmethod =  args$cutmethod

if (is.null(args$rscut)){
  args$rscut <- "0.9"
}
rscut = args$rscut
type = "unsigned"
corType = "pearson"
corFnc = cor
maxPOutliers = ifelse(corType=="pearson",1,0.05)
robustY = ifelse(corType=="pearson",T,F)
## save message
con <- file("WGCNA.AllInOneStep1.log") 
sink(con, append=TRUE) 
sink(con, append=TRUE, type="message") 
## test 
# rawdata ="/Volumes/Samsung_T5/毕业论文/04.WGCNA/01.rawdata/fiber.fpkm.txt"
# samplePerc  = 0.3
# RcCutoff = 1
# datatype =  "FPKM"
# RemainGeneNum = 12000
# method = "rawFPKM"
# cutmethod =  "MAD"
# rscut = 0.9
# datExpr generator -------------------------------------------------------
rawdata = read.table(rawdata,header = T,sep = "\t")
## 01. Filter step1
datExpr01 = getdatExpr(rawdata = rawdata,RcCutoff = RcCutoff,samplePerc = samplePerc,datatype = datatype,method = method)
## 02. Filter step2
GNC = 1-RemainGeneNum/nrow(datExpr01)
datExpr = getdatExpr2(datExpr = datExpr01,GeneNumCut = GNC,cutmethod = cutmethod)
## 03. Get Sample Tree
Stree = getsampleTree(datExpr = datExpr,layout = "circular")
ggsave(Stree$plot,filename = "01.SampleTree.pdf",width = 8,height = 8)
write.tree(Stree$tree,file = "samplecluster.nwk")
## 04. Get Power
beta = getpower(datExpr = datExpr,rscut = rscut)
write.table(beta$sft,file = "01.DetailsofPowerSelection.xls",row.names = F,sep = "\t",quote = F)
print(paste("The Power WGCNA recommand is:",beta$power,sep = " "))
ggsave(beta$plot,filename = "01.SftDetection.pdf",width = 8,height = 6)
save.image(file = "WGCNAoneInAll.Rdata")
# save message ------------------------------------------------------------


sink()
sink(type="message")

cat(readLines("WGCNA.AllInOneStep1.log"), sep="\n")
write.table(cat(readLines("WGCNA.AllInOneStep1.log"), sep="\n"), "log.txt")
